* Clean US PPI, calculate reset prices, perform dispersion regressions

** Updated August 1, 2019 by Lydia Cox to incorporate the BLS PPI micro data. 
	* Current PPI dataset: 2005-2012. 

** Updated August 19, 2019 by Lydia Cox to:
	* Extend PPI dataset back to 1981. 
	* Run new regression specifications. 

** Paths & Folders
global root "/oplcusers2/coxl/PPI"		
global path_logs "$root/logs"
cap mkdir "$path_logs"
cap mkdir "$root/intermediate"
cap mkdir "$root/regressionData"
cap mkdir "$root/savedEstimations"
cap mkdir "$root/graphs"
version 15.1

** Parameters
global exampleData "0"															// Set to zero when applying to real data	
global normaliseGaps "1"														// Demeans and standardise by product
global ageControl "L_xLnPriceAge"												// "" or "L_xLnPriceAge"
 														// Set number of desired bins for bar graphs 

global loadData "1"
	global ld_datasources "usData"										// Options: frData usData
global calcResetPrices "1"														// Also includes competitor prices
	global crp_datasources "usData_dom" 		// Options: frData_dom frData_exp usData_dom usData_exp
	global crp_firstPriceIsReset "0"											// If 1, first observed price is accepted as reset price
global combineIndustries "1"
	global ci_datasources "usData_dom" 		// See crp_datasources
global doRegressions "1"
	global dr_us "1"
		global dru_ebp "1"
		global dru_ff4 "1"
		global dru_bootstrap "0"
global summariseResults "1"
	global sr_14to19 "0"
	global sr_20to25 "1"
	global sr_26 "1"
	global sr_10014 "0"

global horizons "12 24"		// 1 3 12 24

** Define programs
*** dummyGenerator varname, gapcount() middlebinsize() binname() [show graph p() graphbinwidth()]
*** Options, example in parentheses
***		varname: gap variable to process (L_x_Pd)
***		gapcount: number of gaps (7) - middlebinsize: size of middle bin (20) - binname: name of bin variable (gapDummy_res)
***		show: tabulate bin variables - graph: create bar charts - p: price name (Pd)  - graphbinwidth: width of bins in graph (5)
cap program drop dummyGenerator
program define dummyGenerator, rclass
//quietly {
	** Parse options
	syntax varname [aw pw iw fw], gapcount(integer) middlebinsize(real) binname(string) [show] [graph p(string) graphbinwidth(integer 1)]
	local gapvar = "`varlist'"
	local gapCount = "`gapcount'"
	local middleBinSize = "`middlebinsize'"
	local binName = "`binname'"	
	
	** Technicalities
	version 15.1
	tempvar gapDummies
	if "$exampleData" == "1" replace `gapvar' = rnormal(0, 2)
		
	local gapsOneSide = (`gapCount'-1)/2
	local thresholdsOneSide = `gapsOneSide' - 1
	local middleBinSize_neg = 100 - `middleBinSize'
	local middleBinSize_pos = `middleBinSize'
		
	** Cutoffs for middle group
	_pctile `gapvar' [`weight'`exp'] if `gapvar' <= 0, percentiles(`middleBinSize_neg')
	local middleThreshold_neg = r(r1)
	_pctile `gapvar' [`weight'`exp'] if `gapvar' > 0, percentiles(`middleBinSize_pos')
	local middleThreshold_pos = r(r1)

	** Extreme group cutoffs
	_pctile `gapvar' [`weight'`exp'] if `gapvar' < `middleThreshold_neg', nquantiles(`gapsOneSide')
	//mat cutoffsNeg = r(quantiles_used)
	mat cutoffsNeg = J(`thresholdsOneSide',1,.)
	foreach ii of numlist 1/`thresholdsOneSide' {
	mat cutoffsNeg[`ii',1] = r(r`ii')
	} 
	
	_pctile `gapvar' [`weight'`exp'] if `gapvar' > `middleThreshold_pos', nquantiles(`gapsOneSide')
	//mat cutoffsPos = r(quantiles_used)
	mat cutoffsPos = J(`thresholdsOneSide',1,.)
	foreach ii of numlist 1/`thresholdsOneSide' {
		mat cutoffsPos[`ii',1] = r(r`ii')
	} 	

	** Generate variable
	qui gen `gapDummies' = 0 if inrange(`gapvar', `middleThreshold_neg', `middleThreshold_pos')		// Edges inclusive
	qui replace `gapDummies' = 1 if `gapvar' > `middleThreshold_pos' & ~missing(`gapvar')
	qui replace `gapDummies' = -1 if `gapvar' < `middleThreshold_neg'
	qui forvalues i = 1/`thresholdsOneSide' {
		local gapNumber = `i' + 1
		local rowNeg = `thresholdsOneSide' + 1 - `i'
		replace `gapDummies' = `gapNumber'  if `gapvar' > cutoffsPos[`i',1] & ~missing(`gapvar')
		replace `gapDummies' = -`gapNumber' if `gapvar' < cutoffsNeg[`rowNeg',1]
	}
	
	** Rename
	rename `gapDummies' `binName'
	gen `binName'_abs = abs(`binName')
	
	** Show tab
	if "`show'" == "show" {
		noisily tab `binName'
		noisily tab `binName'_abs
	}
	
	** Generate matrix and local of thresholds
	matrix thresholds = cutoffsNeg \ `middleThreshold_neg' \ `middleThreshold_pos' \ cutoffsPos
	foreach i of numlist 1/`=rowsof(thresholds)' {
		local thresholds `thresholds' `=thresholds[`i',1]'
	}
	
	** Make bar charts
	if "`graph'" != "" {
		preserve
		drop if missing(`gapvar')
		tempvar deviation weightvar D_P_up D_P_down
		
		* Round and bound gaps
		if "`graphbinwidth'" != "1" {
			tempvar newgap
			gen double `newgap' = `gapvar'/`graphbinwidth'
			gen double `deviation' = round(`newgap', 0.01)
			replace `deviation' = `deviation'*`graphbinwidth'
		}
		else gen double  `deviation' = round(`gapvar', 0.01)
		sum `deviation', d
		replace `deviation' = r(p99) if `deviation' >= r(p99) & `deviation' != .
		replace `deviation' = r(p1) if `deviation' <= r(p1)
		
		* Calculate statistics
		** Density
		qui sum `deviation', meanonly
		gen `weightvar' = 1
		if "`weight'" != "" replace `weightvar' `exp'
		tab `deviation'
		total `weightvar', over(`deviation')
		mat density = r(table)
		
		** Frequency of price change
		mean D01_`p'_h1 [aw=`weightvar'], over(`deviation')
		mat freq = r(table)
		
		gen `D_P_up' = (DL_`p'_h1 > 0 & DL_`p'_h1 != .)
		mean `D_P_up' [aw=`weightvar'], over(`deviation')
		mat freq_up = r(table)
		
		gen `D_P_down' = (DL_`p'_h1 < 0 & DL_`p'_h1 != .)
		mean `D_P_down' [aw=`weightvar'], over(`deviation')
		mat freq_down = r(table)
		
		** Size of price change, conditional on changeH
		mean DL_`p'_h1 [aw=`weightvar'] if DL_`p'_h1 != 0, over(`deviation')
		mat size = r(table)
		
		mean DL_`p'_h1 [aw=`weightvar'] if (DL_`p'_h1 > 0 & DL_`p'_h1 != .), over(`deviation')
		mat size_up = r(table)
		mean DL_`p'_h1 [aw=`weightvar'] if (DL_`p'_h1 < 0 & DL_`p'_h1 != .), over(`deviation')
		mat size_down = r(table)

		** Get means and SEs into dataset
		local variables = "density freq freq_up freq_down size size_up size_down"
		foreach var of local variables {
			tempfile results`var'
			clear
			mat `var'_t = `var'' 
			xsvmat, from(`var'_t) names(col) fast rownames(`deviation')
			rename * `var'_* 
			drop `var'_t `var'_pvalue `var'_se `var'_df `var'_crit `var'_eform
			rename (`var'_b `var'_ll `var'_ul) (`var'_mean `var'_ll `var'_ul)
			save `results`var''
		} 
		
		clear
		gen `deviation' = ""
		foreach var of local variables {
			merge 1:1 `deviation' using `results`var'', nogen
		}
		
		** Graph
		drop if real(`deviation') == .
		destring `deviation', replace
		label var `deviation' "Gap"
		foreach var of local variables {
			twoway (bar `var'_mean `deviation', barwidth(`=`graphbinwidth'/100') fcolor(black%50) lcolor(white))(rcap `var'_ll `var'_ul `deviation', lcolor(black%80)), xline(`thresholds', lcolor(black%25) lpattern(dash)) title(`gapvar' - `var')
			graph export "$root/graphs/`gapvar'_`var'.png", replace width(1600)
			graph export "$root/graphs/`gapvar'_`var'.pdf", replace
			graph save "$root/graphs/`gapvar'_`var'.gph", replace
		}
	}
	
	** Return matrix and local of thresholds
	noisily matrix list thresholds
	return matrix thresholds = thresholds
//}
end

*** dummyRegs, depvar() dummyvar() shock() [shock_interaction(varname)] controls() absorb() cluster() [regnumber()]
*** Note: will loop over up/down/abs and horizons
*** Options, example in parentheses
***		depvar: dependent variable without h, with @ where up/down is (D_Pd_@)
***		dummyvar: dummy variable (gapDummy_res) - shock: shock variable (ebpnew) - shock_interaction: optional interaction version of shock (ebpnew_residualP)
***		controls: ... (L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa) - absorb: FE (month12 productID) - cluster: ... (month industryID)
***		regnumber: regression number to save estimates, optional (18ff4)
cap program drop dummyRegs
program define dummyRegs
	* Parse options
	syntax [if] [in] [aw pw iw fw], depvar(string) dummyvar(varname) shock(varname) [shock_interaction(varname)] controls(varlist) absorb(varlist) cluster(varlist) [regnumber(string)]
	if "`shock_interaction'" == "" local shock_interaction = "`shock'"
	
	// 
	tempvar binString
	tostring `dummyvar', gen(`binString')
	replace `binString' = subinstr(`binString', "-", "M", .)							// Replace - by M
	replace `binString' = "P" + `binString' if (substr(`binString', 1, 1) != "M" & `binString' != ".")	// Add P for positive values
	replace `binString' = "" if `binString' == "." 
	
	* Prepare dummies and interactions
	** Regular
	levelsof `binString', local(bins)
	foreach bin of local bins {
		gen `dummyvar'_`bin' = (`binString' == "`bin'") if `dummyvar' != .		// Split categorical into dummy variables
		gen `shock_interaction'_`bin' = `dummyvar'_`bin' * `shock_interaction'						// Interact dummies with shock
	}
	
	** Absolute values
	tempvar binString_abs
	levelsof `dummyvar'_abs, local(bins)
	foreach bin of local bins {
		gen `dummyvar'_abs_`bin' = (`dummyvar'_abs == `bin') if `dummyvar'_abs != .
		gen `shock_interaction'_abs_`bin' = `dummyvar'_abs_`bin' * `shock_interaction'
	}
	
	local depvar_up = subinstr("`depvar'", "@", "up", 1)
	local depvar_down = subinstr("`depvar'", "@", "down", 1)
	local depvar_abs = subinstr("`depvar'", "_@", "", 1)
	
	* Run regressions
	drop `dummyvar'_P0 `shock_interaction'_P0 `shock_interaction'_abs_0
	*drop `dummyvar'_M3 `shock_interaction'_M3		// To compare to Stata way
	
	foreach h of global horizons {
		/* up */ 	reghdfe `depvar_up'_h`h' `dummyvar'_M* `dummyvar'_P* `shock_interaction'_M* `shock_interaction'_P* `shock' `controls' [`weight'`exp'] `if' `in', absorb(`absorb') cluster(`cluster')
		noisily di "Above results: `e(cmdline)'"
		estimates save "$root/savedEstimations/`regnumber'up_h`h'", replace
		/* down */ 	reghdfe `depvar_down'_h`h' `dummyvar'_M* `dummyvar'_P* `shock_interaction'_M* `shock_interaction'_P* `shock' `controls' [`weight'`exp'] `if' `in', absorb(`absorb') cluster(`cluster')
		noisily di "Above results: `e(cmdline)'"
		estimates save "$root/savedEstimations/`regnumber'down_h`h'", replace
		/*/* abs */ 	reghdfe `depvar_abs'_h`h' `dummyvar'_M* `dummyvar'_P* `shock_interaction'_abs_* `shock' `controls' [`weight'`exp'] `if' `in', absorb(`absorb') cluster(`cluster')
		noisily di "Above results: `e(cmdline)'"
		estimates save "$root/savedEstimations/`regnumber'abs_h`h'", replace
		*/
	}
	
	drop `dummyvar'_M* `dummyvar'_P* `shock_interaction'_M* `shock_interaction'_P* `dummyvar'_abs_* `shock_interaction'_abs_*
end

*** baseRegs: See dummyRegs, but with two added options: gapvar(L_x_P) and absgapvar(absGap_res)
*** Note: Can run both base regressions and time fixed effects ones (simply adjust the options)
cap program drop baseRegs
program define baseRegs
	* Parse options
	syntax [if] [in] [aw pw iw fw], depvar(string) gapvar(varname) absgapvar(varname) [shock(varname) shock_interaction(varname)] controls(varlist) absorb(varlist) cluster(varlist) [regnumber(string)]
	if "`shock_interaction'" == "" local shock_interaction = "`shock'"
	
	* Depvars
	local depvar_up = subinstr("`depvar'", "@", "up", 1)
	local depvar_down = subinstr("`depvar'", "@", "down", 1)
	local depvar_abs = subinstr("`depvar'", "_@", "", 1)
	
	* Run Regressions
	foreach h of global horizons {
		/* up */ 	reghdfe `depvar_up'_h`h' 	`gapvar' c.`gapvar'#c.`shock_interaction' `shock' `controls' [`weight'`exp'] `if' `in', absorb(`absorb') cluster(`cluster')
		noisily di "Above results: `e(cmdline)'"
		estimates save "$root/savedEstimations/`regnumber'up_h`h'", replace
		/* down */ 	reghdfe `depvar_down'_h`h' 	`gapvar' c.`gapvar'#c.`shock_interaction' `shock' `controls' [`weight'`exp'] `if' `in', absorb(`absorb') cluster(`cluster')
		noisily di "Above results: `e(cmdline)'"
		estimates save "$root/savedEstimations/`regnumber'down_h`h'", replace
		/*/* abs */	reghdfe `depvar_abs'_h`h' 	`gapvar' c.`absgapvar'#c.`shock_interaction' `shock' `controls' [`weight'`exp'] `if' `in', absorb(`absorb') cluster(`cluster')
		noisily di "Above results: `e(cmdline)'"
		estimates save "$root/savedEstimations/`regnumber'abs_h`h'", replace
		*/
	}
end

*** allRegs, See dummyRegs and baseRegs for depvar/shock/shock°interaction/controls/absorb/cluster
***		gapvars: reset and competitor gap vars (L_x_Pd L_xCompDem_Pd)
***		absgapvars/dummyvars: same system as gapvars
***		controls_tfe/absorb_tfe: TFE versions
***		regnumbers: separate with whitespace (14 16 18)
cap program drop allRegs
program define allRegs
	* Parse options
	syntax [if] [in] [aw pw iw fw], depvar(string) gapvars(varlist) absgapvars(varlist) dummyvars(varlist) [shock(varname) shock_interaction(varname)] controls(varlist) controls_tfe(varlist) absorb(varlist) absorb_tfe(varlist) cluster(varlist) [regnumbers(string)]
	if "`shock_interaction'" == "" local shock_interaction = "`shock'"
	local gapvarR : word 1 of `gapvars'
	local gapvarC : word 2 of `gapvars'
	local absgapvarR : word 1 of `absgapvars'
	local absgapvarC : word 2 of `absgapvars'
	local dummyvarR : word 1 of `dummyvars'
	local dummyvarC : word 2 of `dummyvars'
	
	local regnumber_base 	: word 1 of `regnumbers'
	local regnumber_tfe 	: word 2 of `regnumbers'
	local regnumber_dummy 	: word 3 of `regnumbers'
	
	* Run regressions
	** Base
	noisily di "Base regs are numbered `regnumber_base' and 10`regnumber_base'"
	baseRegs `if' `in' [`weight'`exp'], depvar(`depvar') gapvar(`gapvarR') absgapvar(`absgapvarR') 	shock(`shock') shock_interaction(`shock_interaction') controls(`controls') absorb(`absorb') cluster(`cluster') regnumber(`regnumber_base')
	baseRegs `if' `in' [`weight'`exp'], depvar(`depvar') gapvar(`gapvarC') absgapvar(`absgapvarC') 	shock(`shock') shock_interaction(`shock_interaction') controls(`controls') absorb(`absorb') cluster(`cluster') regnumber(10`regnumber_base')
	
	** TFE
	noisily di "TFE regs are numbered `regnumber_tfe' and 10`regnumber_tfe'"
	baseRegs `if' `in' [`weight'`exp'], depvar(`depvar') gapvar(`gapvarR') absgapvar(`absgapvarR') 	shock_interaction(`shock_interaction') controls(`controls_tfe') absorb(`absorb_tfe') cluster(`cluster') regnumber(`regnumber_tfe')
	baseRegs `if' `in' [`weight'`exp'], depvar(`depvar') gapvar(`gapvarC') absgapvar(`absgapvarC') 	shock_interaction(`shock_interaction') controls(`controls_tfe') absorb(`absorb_tfe') cluster(`cluster') regnumber(10`regnumber_tfe')
	
	** Dummy
	noisily di "Dummy regs are numbered `regnumber_dummy' and 10`regnumber_dummy'"
	dummyRegs `if' `in' [`weight'`exp'], depvar(`depvar') dummyvar(`dummyvarR') shock(`shock') shock_interaction(`shock_interaction') controls(`controls') absorb(`absorb') cluster(`cluster') regnumber(`regnumber_dummy')
	dummyRegs `if' `in' [`weight'`exp'], depvar(`depvar') dummyvar(`dummyvarC') shock(`shock') shock_interaction(`shock_interaction') controls(`controls') absorb(`absorb') cluster(`cluster') regnumber(10`regnumber_dummy')
end
				
** Code
if "$loadData" == "1" {
    * Start log
    cap log close loadData
    log using "$path_logs/b1_loadData.log", replace text name(loadData)
	
	if strpos("$ld_datasources", "frData") {
		* Start log
		cap log close ld_datasource
		log using "$path_logs/b1_ld_datasource_${ld_datasource}.log", replace text name(ld_datasource)
		
		* Load data
		use "$root/input/export_data_new", clear

		* Clean up a bit
		** Variables
		rename (price_D price_X unique_ID firmid item_wt month) (price_dom price_exp productID firmID weight_item month_nr)
		gen month = ym(year, month_nr)
		
		order ind_4d firmID productID year month price_exp price_dom weight_item
		format month %tm
		drop month_nr
		
		sort ind_4d firmID productID year month 
		
		** Observations
		duplicates drop ind_4d firmID productID year month, force
		
		** Missings
		*** Prep
		xtset productID month
		tsfill
		
		*** Export price missings
		gen gap_exp = 1 if ~missing(price_exp[_n-1], price_exp[_n+1]) & missing(price_exp)
		replace price_exp = L.price_exp if gap_exp == 1
		replace ind_4d = ind_4d[_n-1] if gap_exp == 1
		replace firmID = firmID[_n-1] if gap_exp == 1
		drop gap_exp
		
		*** Domestic price missings
		gen gap_dom = 1 if ~missing(price_dom[_n-1], price_dom[_n+1]) & missing(price_dom)
		replace price_dom = L.price_dom if gap_dom == 1
		replace ind_4d = ind_4d[_n-1] if gap_dom == 1
		replace firmID = firmID[_n-1] if gap_dom == 1
		drop gap_dom
		
		* Numerify industries
		gen industryID = real(ind_4d)
		order industryID, after(ind_4d)
		drop ind_4d

		
		* Save data
		tempfile allfrData
		save `allfrData'
		
		** Export prices
		*** Reduce to export prices
		drop price_dom
		rename price_exp price
		save "$root/intermediate/frData_exp_1_clean", replace
		
		*** Split into categories
		levelsof industryID
		foreach ind in `r(levels)' {
			use "$root/intermediate/frData_exp_1_clean", replace
			keep if industryID == `ind'
			save "$root/intermediate/frData_exp_1_clean_`ind'", replace
		}
		
		** Domestic prices
		*** Reduce to domestic prices
		use `allfrData', clear
		drop price_exp
		rename price_dom price
		save "$root/intermediate/frData_dom_1_clean", replace
		
		*** Split into categories
		levelsof industryID
		foreach ind in `r(levels)' {
			use "$root/intermediate/frData_dom_1_clean", replace
			keep if industryID == `ind'
			save "$root/intermediate/frData_dom_1_clean_`ind'", replace
		}

		** Restore full data
		use `allfrData', clear
		
		* Close log
		log close ld_datasource
	}
	
	if strpos("$ld_datasources", "usData") {
		* Start log
		cap log close ld_datasource
		log using "$path_logs/b1_ld_datasource_${ld_datasource}.log", replace text name(ld_datasource)
		
		* PPI
		** Load data
		if "$exampleData" == "1" use "$root/Input/fakeppi_v1", clear
		else use "$root/input/usData_PPI_81to12.dta", clear

		** Clean up a bit
		*** Variables
		rename (prc item_code firmid item_wgt rel_wgt shipments naics_item) (price productID firmID weight_withinFirm weight_betweenFirm weight_item ind_6d)
		drop dlprc
		drop if ind_6d == . 
		drop weight_*Firm
		
		/* sas date is daily, convert this to monthly */ 
		gen month = mofd(date)
		format month %tm 
		gen year = year(dofm(month))
		drop date
		
		sort ind_6d firmID productID year month 
		
		** Observations
		duplicates drop ind_6d firmID productID year month, force

		** Missings
		*** Prep
		
		//productID Is a string, and can't xtset strings. Too many values to encode, so we create a numeric code
		//and save string code as productID_string
		egen productIDcode = group(productID)
		rename productID productID_string
		rename productIDcode productID
		xtset productID month
		tsfill
		
		*** Export price missings
		gen gap = 1 if ~missing(price[_n-1], price[_n+1]) & missing(price)
		replace price = L.price if gap == 1
		replace ind_6d = ind_6d[_n-1] if gap == 1
		replace firmID = firmID[_n-1] if gap == 1
		drop gap
		
		** Generate aggregated categories
		tostring(ind_6d), replace 
		gen ind_4d = substr(ind_6d, 1, 4)		// Change here for different industry aggregation (1/2) #indClass

		* Numerify industries
		gen industryID = real(ind_4d)			// (not necessary to change, but could be nice to use ind_3d instead of ind_4d to avoid confusion)
		order industryID, after(ind_4d)			// ...
		drop ind_4d ind_6d						// ...

		
		** Save data
		save "$root/intermediate/usData_dom_1_clean", replace
		
		*** Split into categories
		preserve
		levelsof industryID
		foreach ind in `r(levels)' {
			use "$root/intermediate/usData_dom_1_clean", replace
			keep if industryID == `ind'
			save "$root/intermediate/usData_dom_1_clean_`ind'", replace
		}
		restore
		
		* IPP
		** Load data
		if "$exampleData" == "1" {
			use "$root/Input/fakeipp_v1", clear
			gen weight_item = 1
		}
		else use "$root/input/usData_IPP", clear
		gen weight_item = 1 
		
		** Clean up a bit
		*** Variables
		rename (item_code rts_id hts6 net_price) (productID firmID ind_6d price) // 
		drop dlprc
		
		gen month = mofd(mydate)
		format month %tm
		gen year = year(dofm(month))
		drop mydate
		
		sort ind_6d firmID productID year month 
		
		** Observations
		duplicates drop ind_6d firmID productID year month, force

		** Missings
		*** Prep
		egen productIDcode = group(productID)
		rename productID productID_string
		rename productIDcode productID
		xtset productID month, monthly
		tsfill
		
		*** Export price missings
		gen gap = 1 if ~missing(price[_n-1], price[_n+1]) & missing(price)
		replace price = L.price if gap == 1
		replace ind_6d = ind_6d[_n-1] if gap == 1
		replace firmID = firmID[_n-1] if gap == 1
		drop gap

		** Generate aggregated categories
		tostring(ind_6d), replace 
		gen ind_4d = substr(ind_6d, 1, 4) 		// Change here for different industry aggregation (2/2) #indClass
		
		* Numerify industries
		gen industryID = real(ind_4d)			// (not necessary to change, but could be nice to use ind_3d instead of ind_4d to avoid confusion)
		order industryID, after(ind_4d)			// ...
		drop ind_4d ind_6d						// ...
	
		** Save data
		save "$root/intermediate/usData_exp_1_clean", replace
		
		*** Split into categories
		preserve
		levelsof industryID
		foreach ind in `r(levels)' {
			use "$root/intermediate/usData_exp_1_clean", replace
			keep if industryID == `ind'
			save "$root/intermediate/usData_exp_1_clean_`ind'", replace
		}
		restore
		
		* Close log
		log close ld_datasource
	}
    
    * Close log
    log close loadData
}

if "$calcResetPrices" == "1" {
    * Start log
    cap log close calcResetPrices
    log using "$path_logs/b1_calcResetPrices.log", replace text name(calcResetPrices)
	
	* Loop over datasources
	foreach datasource of global crp_datasources {
		di as result _newline _newline "`datasource'" _newline
		
		* Loop over industries
		local allFiles : dir "$root/intermediate/" files "`datasource'_1_clean_*.dta", respectcase
		qui foreach file of local allFiles {
			* Load data 
			local industry = subinstr("`file'", "`datasource'_1_clean_", "", .)
			local industry = subinstr("`industry'", ".dta", "", .)
			use "$root/intermediate/`datasource'_1_clean_`industry'", replace
			
			* Reset Price Calculation
			noisily di "`industry'"
			
			** Prep
			*** Generate log prices
			gen P = log(price)
			drop price

			*** Identify price changes (relative to previous observation)
			bysort productID (month): gen byte D_P = (P != P[_n-1]) if _n != 1 & missing(P, P[_n-1]) == 0

			*** Reset price equals current price if the price changed 
			gen P_reset = P if D_P == 1
			if "$crp_firstPriceIsReset" == "1" bysort productID (month): replace P_reset = P if _n == 1

			*** Weight used to calculate reset price inflation (nonmissing for price changers only)
			bysort month (productID): egen tot_weight_r = total(weight_item) if D_P == 1
			gen P_weight_reset = weight_item/tot_weight_r if D_P == 1
			drop tot_weight_r

			** Reset price & inflation calculation
			matrix drop _all
			sort productID month
			sum month
			matrix DL_P_reset = (r(min),.,.)
			foreach month of numlist `=r(min) + 1'/`r(max)' {
				* Show progress
				noisily di %4s = "`month' " _continue
				if `month'/30 == trunc(`month'/30) noisily di " |"

				* Calculate product level reset price inflation
				gen DL_P_reset_productID = D.P_reset	if month == `month'
				
				* Calculate weighted average category level reset price inflation
				sum DL_P_reset_productID [aweight=P_weight_reset], meanonly
				local DL_P_reset_`month' = r(sum)
				matrix DL_P_reset = (DL_P_reset \ (`month', r(sum), r(N)))

				* Update reset prices of productIDs that didn't change their price
				replace P_reset = L.P_reset + `DL_P_reset_`month''	if D_P != 1	& month == `month' & ~missing(L.P_reset)

				* Clean up and restart
				drop DL_P_reset_productID
			}
			noisily di _newline
			
			** Clean up
			drop P_weight_reset
			
			** Obtain reset price inflations
			svmat DL_P_reset
			rename DL_P_reset1 DL_P_reset_month
			rename DL_P_reset2 DL_P_reset
			rename DL_P_reset3 D_P_count
			format DL_P_reset_month %tm
			
			* Calculate competitor prices
			bysort industryID month: egen total_PRF = total(P)
			bysort industryID month: egen count_PRF = count(P)
			gen meanOthers_P = (total_P - P)/(count_P - 1)
			drop total_P
			
			** Calculate deviation from competitor prices
			gen gapComp_P = P - meanOthers_P
			drop meanOthers_P
			
			** Demean by firm average (store in IRI)
			bysort firmID: egen mean_gapComp_P = mean(gapComp_P)
			gen gapCompDemeanedFirmInd_P = gapComp_P - mean_gapComp_P
			drop mean_gapComp_P
			
			* Lag
			xtset productID month, noquery
			sort productID month
			
			gen L_xComp_P = L.gapComp_P
			gen L_xCompDem_P = L.gapCompDemeanedFirmInd_P
			
			label variable L_xComp_P "Lagged gap own price vs average competitor price"
			label variable L_xCompDem_P "Lagged gap own price vs average competitor price. Demeaned by firm X industry."

			* Save
			save "$root/intermediate/`datasource'_2_wResetPrices_`industry'", replace
		}
	}
	
    * Close log
    log close calcResetPrices
}
******************** COMBINE INDUSTRIES ************************
if "$combineIndustries" == "1" {
    * Start log
    cap log close combineIndustries
    log using "$path_logs/b1_combineIndustries.log", replace text name(combineIndustries)
	
	* Loop over data sources
	foreach datasource of global ci_datasources {
		di as result _newline _newline "`datasource'" _newline
		
		* Append industries
		clear
		local allFiles : dir "$root/intermediate/" files "`datasource'_2_wResetPrices_*.dta", respectcase
		qui foreach file of local allFiles {
			append using "$root/intermediate/`file'", nolabel
		}
		
		* Calculate horizon & split price changes
		** Split price increases and decreases, h1
		sort productID month
		rename D_P D_P_h1
		gen DL_P_h1 = P - L.P
		gen D01_P_h1 = D_P_h1
		replace D_P_h1 = -D_P_h1 if DL_P_h1 < 0
		
		gen D_P_up_h1 = 1 		if D_P_h1 == 1
		replace D_P_up_h1 = 0 	if inlist(D_P_h1, 0, -1)
		gen D_P_down_h1 = 1 	if D_P_h1 == -1
		replace D_P_down_h1 = 0 if inlist(D_P_h1, 0, 1)
		
		** Calculate 12 month price changes
		gen DL_P_h12 = F11.P - L.P
		gen D_P_h12 = sign(DL_P_h12)
		gen D01_P_h12 = abs(D_P_h12)
		
		gen D_P_up_h12 = 1 		if D_P_h12 == 1
		replace D_P_up_h12 = 0 	if inlist(D_P_h12, 0, -1)
		gen D_P_down_h12 = 1 	if D_P_h12 == -1
		replace D_P_down_h12 = 0 if inlist(D_P_h12, 0, 1)
		
		** Calculate 12 month price changes
		gen DL_P_h24 = F23.P - L.P
		gen D_P_h24 = sign(DL_P_h24)
		gen D01_P_h24 = abs(D_P_h24)
		
		gen D_P_up_h24 = 1 		if D_P_h24 == 1
		replace D_P_up_h24 = 0 	if inlist(D_P_h24, 0, -1)
		gen D_P_down_h24 = 1 	if D_P_h24 == -1
		replace D_P_down_h24 = 0 if inlist(D_P_h24, 0, 1)
		
		** Calculate 3 month price changes
		gen DL_P_h3 = F2.P - L.P
		gen D_P_h3 = sign(DL_P_h3)
		gen D01_P_h3 = abs(D_P_h3)
		
		gen D_P_up_h3 = 1 		if D_P_h3 == 1
		replace D_P_up_h3 = 0 	if inlist(D_P_h3, 0, -1)
		gen D_P_down_h3 = 1 	if D_P_h3 == -1
		replace D_P_down_h3 = 0 if inlist(D_P_h3, 0, 1)
		
		** Lagged price gap
		gen double x_P = P - P_reset
		gen double L_x_P = L.x_P
		
		* Merge in shocks
		/* FF4 */ merge m:1 month using "$root/input/Karadi2017.dta", keepusing(FF4_alt1_norm) keep(match master) nogen
		/* ECB */ merge m:1 month using "$root/input/Karadi2019_ECB.dta", keepusing(mpshockpoorman) keep(match master) nogen
		sort productID month
		
		** Generate gap dummies
		*** Reset prices
		sum L_x_P if L_x_P < 0, d
		local median_neg = r(p50)												// In time we might want to hardcode this to avoid issues
		local p1_neg = r(p99)												// In time we might want to hardcode this to avoid issues
		sum L_x_P if L_x_P > 0, d
		local median_pos = r(p50)												// In time we might want to hardcode this to avoid issues
		local p1_pos = r(p1)												// In time we might want to hardcode this to avoid issues
		gen 	P_Lgap = 1 	if inrange(L_x_P, -9999, `median_neg')
		replace P_Lgap = 2 	if inrange(L_x_P, `median_neg', 0)
		replace P_Lgap = 3 	if L_x_P == 0
		replace P_Lgap = 4 	if inrange(L_x_P, 0, `median_pos')
		replace P_Lgap = 5 	if inrange(L_x_P, `median_pos', 9999)
		
		*** Competitor prices
		sum L_xCompDem_P if L_xCompDem_P < 0, d
		local median_neg = r(p50)
		sum L_xCompDem_P if L_xCompDem_P > 0, d
		local median_pos = r(p50)
		gen 	P_LgapComp = 1 	if inrange(L_xCompDem_P, -9999, `median_neg')
		replace P_LgapComp = 2 	if inrange(L_xCompDem_P, `median_neg', 0)
		replace P_LgapComp = 3 	if L_xCompDem_P == 0
		replace P_LgapComp = 4 	if inrange(L_xCompDem_P, 0, `median_pos')
		replace P_LgapComp = 5 	if inrange(L_xCompDem_P, `median_pos', 9999)
		
		label define P_Lgap 1 "Negative, Large" 2 "Negative, Small" 3 "Zero" 4 "Positive, Small" 5 "Positive, Large"
		label values P_Lgap P_LgapComp P_Lgap
		tab P_Lgap
		tab P_LgapComp

		* Save
		*drop if missing(L_x_P)
		save "$root/regressionData/`datasource'_3_prepped", replace
	}
    
    * Close log
    log close combineIndustries
}

if "$doRegressions" == "1" {
    * Start log
    cap log close doRegressions
    log using "$path_logs/b1_doRegressions.log", replace text name(doRegressions)
	
	if "$dr_us" == "1" {
		* Start log
		cap log close dr_us
		log using "$path_logs/b1_dr_us.log", replace text name(dr_us)
		
		* Load domestic data
		use "$root/regressionData/usData_dom_3_prepped", replace
		rename *_P_* *_Pd_*
		rename *_P *_Pd
		rename P_* Pd_*
		rename P Pd
		gen month12 = month(dofm(month))
		
		* Merge in extra variables
		merge m:1 month using "$root/input/US_variables.dta", nogen keep(match master using)		
		drop L7_* L8_* L9_* L10_* L11_* L12_*
		sort productID month
		
		* Normalise
		// NEW version of normalization 
		/*if "$normaliseGaps" == "1" {
		
			tempvar pricechangeSize
			
			gen `pricechangeSize' = DL_Pd_h1 if abs(DL_Pd_h1) >= 0.01
			bysort productID: center L_x_Pd L_xCompDem_Pd, inplace 
			bysort productID: egen DL_Pd_h1_mean = mean(abs(`pricechangeSize'))
			foreach var of varlist L_x_Pd L_xCompDem_Pd {
				replace `var' = `var'/DL_Pd_h1_mean if DL_Pd_h1_mean != 0
				bysort productID: egen temp = sd(abs(`var'))
				replace `var' = `var'/temp if DL_Pd_h1_mean == 0
				drop temp
			}
				drop DL_Pd_h1_mean
		}
		*/
		
		// OLD version of normalization 
		if "$normaliseGaps" == "1" {
			bysort productID: center L_x_Pd L_xCompDem_Pd, inplace
			bysort productID: egen DL_Pd_h1_mean = mean(abs(DL_Pd_h1))
			foreach var of varlist L_x_Pd L_xCompDem_Pd {
				replace `var' = `var'/DL_Pd_h1_mean if DL_Pd_h1_mean != 0
				
				bysort productID: egen temp = sd(abs(`var'))
				replace `var' = `var'/temp if DL_Pd_h1_mean == 0
				drop temp
			}
			drop DL_Pd_h1_mean
		}
		
		* Identify high and low frequency industries
		bysort industryID: egen frequency_industryMean = mean(D01_Pd_h1)
		egen tag = tag(industryID)
		sum frequency_industryMean if tag == 1, d
		gen industry_highFreq = (frequency_industryMean >= r(p50))
		tab industryID industry_highFreq
		drop frequency_industryMean tag
		
		* Generate price age		// Requires no missing values (gaps are okay, but need to be pre-dropped)
		** Identify price change moments
		gsort productID month
		gen DP_month = month if D01_Pd_h1 == 1
		by productID (month): replace DP_month = DP_month[_n-1] if ~missing(DP_month[_n-1]) & missing(DP_month)
		
		** Months since price change
		gen priceAge = month - DP_month + 1
		gen L_xPriceAge = L.priceAge
		gen L_xLnPriceAge = log(L_xPriceAge)
		
		* Split gaps into groups 
		
		** Calculate bin size 
		local nBins = 31 
		sum L_x_Pd, d
		local binwidth = int(round(100*(r(p99)-r(p1))/`nBins'))
		dummyGenerator L_x_Pd, gapcount(7) middlebinsize(10) binname(gapDummy_res) show graph p(Pd) graphbinwidth(`binwidth')
		
		** Calculate bin size 
		sum L_xCompDem_Pd, d
		local binwidth = int(round(100*(r(p99)-r(p1))/`nBins'))
		dummyGenerator L_xCompDem_Pd, gapcount(7) middlebinsize(10) binname(gapDummy_comp) show graph p(Pd) graphbinwidth(`binwidth')
		
		
		* Alternative gap measures
		gen positiveGap = (L_x_P >= 0) & L_x_P != .
		gen absGap_res = abs(L_x_P)
		gen absGap_comp = abs(L_xCompDem_P)
			

		* Run regressions
		if "$dru_ebp" == "1" {
			* Start log
			cap log close dru_ebp
			log using "$path_logs/b1_dru_ebp_dom.log", replace text name(dru_ebp)
			
			* Main
			allRegs, depvar(D_Pd_@) gapvars(L_x_Pd L_xCompDem_Pd) absgapvars(absGap_res absGap_comp) dummyvars(gapDummy_res gapDummy_comp) shock(ebpnew) shock_interaction(ebpnew_residualP) controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}) controls_tfe(${ageControl}) absorb_tfe(month productID) absorb(month12 productID) cluster(month industryID) regnumbers(14 16 18)

			* Robustness
			** Pre 2007
			baseRegs if year <= 2007, depvar(D_Pd_@) gapvar(L_x_P) absgapvar(absGap_res) 			shock(ebpnew) shock_interaction(ebpnew_residualP) controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(20)
			baseRegs if year <= 2007, depvar(D_Pd_@) gapvar(L_xCompDem_P) absgapvar(absGap_comp) 	shock(ebpnew) shock_interaction(ebpnew_residualP) controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(1020)
			
			** Low and high frequency
			baseRegs if industry_highFreq == 0, depvar(D_Pd_@) gapvar(L_x_P) absgapvar(absGap_res) 			shock(ebpnew) shock_interaction(ebpnew_residualP) controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(22)
			baseRegs if industry_highFreq == 0, depvar(D_Pd_@) gapvar(L_xCompDem_P) absgapvar(absGap_comp) 	shock(ebpnew) shock_interaction(ebpnew_residualP) controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(1022)
			baseRegs if industry_highFreq == 1, depvar(D_Pd_@) gapvar(L_x_P) absgapvar(absGap_res) 			shock(ebpnew) shock_interaction(ebpnew_residualP) controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(24)
			baseRegs if industry_highFreq == 1, depvar(D_Pd_@) gapvar(L_xCompDem_P) absgapvar(absGap_comp) 	shock(ebpnew) shock_interaction(ebpnew_residualP) controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(1024)
			
			** With weights	//#weightregression1
			baseRegs [aw=weight_item], depvar(D_Pd_@) gapvar(L_x_P) absgapvar(absGap_res) 			shock(ebpnew) shock_interaction(ebpnew_residualP) controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(26)
			baseRegs [aw=weight_item], depvar(D_Pd_@) gapvar(L_xCompDem_P) absgapvar(absGap_comp) 	shock(ebpnew) shock_interaction(ebpnew_residualP) controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(1026)
			
			* Close log
			log close dru_ebp
		}
		
		if "$dru_ff4" == "1" {
			* Start log
			cap log close dru_ff4
			log using "$path_logs/b1_dru_ff4_dom.log", replace text name(dru_ff4)
			label var FF4_alt1_norm "ff4"
			
			* Main
			allRegs, depvar(D_Pd_@) gapvars(L_x_Pd L_xCompDem_Pd) absgapvars(absGap_res absGap_comp) dummyvars(gapDummy_res gapDummy_comp) shock(FF4_alt1_norm) controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew ${ageControl}) controls_tfe(${ageControl}) absorb_tfe(month productID) absorb(month12 productID) cluster(month industryID) regnumbers(14ff4_ 16ff4_ 18ff4_)

			* Robustness
			** Pre 2007
			baseRegs if year <= 2007, depvar(D_Pd_@) gapvar(L_x_P) absgapvar(absGap_res) 			shock(FF4_alt1_norm)  controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(20ff4_)
			baseRegs if year <= 2007, depvar(D_Pd_@) gapvar(L_xCompDem_P) absgapvar(absGap_comp) 	shock(FF4_alt1_norm)  controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(1020ff4_)
			
			** Low and high frequency
			baseRegs if industry_highFreq == 0, depvar(D_Pd_@) gapvar(L_x_P) absgapvar(absGap_res) 			shock(FF4_alt1_norm)  controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(22ff4_)
			baseRegs if industry_highFreq == 0, depvar(D_Pd_@) gapvar(L_xCompDem_P) absgapvar(absGap_comp) 	shock(FF4_alt1_norm)  controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(1022ff4_)
			baseRegs if industry_highFreq == 1, depvar(D_Pd_@) gapvar(L_x_P) absgapvar(absGap_res) 			shock(FF4_alt1_norm)  controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(24ff4_)
			baseRegs if industry_highFreq == 1, depvar(D_Pd_@) gapvar(L_xCompDem_P) absgapvar(absGap_comp) 	shock(FF4_alt1_norm)  controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(1024ff4_)
			
			** With weights 
			baseRegs [aw=weight_item], depvar(D_Pd_@) gapvar(L_x_P) absgapvar(absGap_res) 			shock(FF4_alt1_norm)  controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(26ff4_)
			baseRegs [aw=weight_item], depvar(D_Pd_@) gapvar(L_xCompDem_P) absgapvar(absGap_comp) 	shock(FF4_alt1_norm)  controls(L*_lip L*_lppi_sa L*_gs1 L*_ebpnew ${ageControl}) absorb(month12 productID) cluster(month industryID) regnumber(1026ff4_)

						
			* Close log
			log close dru_ff4
		}

		if "$dru_bootstrap" == "1" {
			* Start log
			cap log close dru_bootstrap
			log using "$path_logs/b1_dru_bootstrap.log", replace text name(dru_bootstrap)
			
			* Predict ebpnew_residual
			egen tag = tag(month)
			reg ebpnew L*_gs1 L*_lip L*_lppi_sa L*_ebpnew if tag == 1, r
			predict ebpnew_residualP_bs, r
			compare ebpnew_residualP*
			drop tag ebpnew_residualP_bs
			
			* Do one regression
			reghdfe D_Pd_up_h12 L_x_P c.L_x_P#c.ebpnew_residual ebpnew L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}, absorb(month12 productID) cluster(month industryID)
			estimates save "$root/savedEstimations/10014up_LPM_FEidnr12_h12_Pd", replace
			reghdfe D_Pd_down_h12 L_x_P c.L_x_P#c.ebpnew_residual ebpnew L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}, absorb(month12 productID) cluster(month industryID)
			estimates save "$root/savedEstimations/10014down_LPM_FEidnr12_h12_Pd", replace
			
			* Program define
			xtset, clear
			cap program drop bootstrap_14u
			program define bootstrap_14u
				cap drop tag ebpnew_residualP_bs
				egen tag = tag(month_bs)
				reg ebpnew L*_gs1 L*_lip L*_lppi_sa L*_ebpnew if tag == 1, r
				predict ebpnew_residualP_bs, r
				
				reghdfe D_Pd_up_h12 L_x_P c.L_x_P#c.ebpnew_residualP_bs ebpnew L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}, absorb(month12 productID) cluster(month_bs industryID)
			end
			
			cap program drop bootstrap_14d
			program define bootstrap_14d
				cap drop tag ebpnew_residualP_bs
				egen tag = tag(month_bs)
				reg ebpnew L*_gs1 L*_lip L*_lppi_sa L*_ebpnew if tag == 1, r
				predict ebpnew_residualP_bs, r
				
				reghdfe D_Pd_down_h12 L_x_P c.L_x_P#c.ebpnew_residualP_bs ebpnew L*_lip L*_lppi_sa L*_gs1 L*_ebpnew gs1 lip lppi_sa ${ageControl}, absorb(month12 productID) cluster(month_bs industryID)
			end
			
			bootstrap, reps(100) cluster(month) idcluster(month_bs): bootstrap_14u
			estimates save "$root/savedEstimations/20014up_LPM_FEidnr12_h12_Pd", replace
			bootstrap, reps(100) cluster(month) idcluster(month_bs): bootstrap_14d
			estimates save "$root/savedEstimations/20014down_LPM_FEidnr12_h12_Pd", replace
			
			* Close log
			log close dru_bootstrap
		}

		* Close log
		log close dr_us
	}

    * Close log
   log close doRegressions
}

if "$summariseResults" == "1" {
    * Start log
    cap log close summariseResults
    log using "$path_logs/b1_summariseResults.log", replace text name(summariseResults)
	set linesize 250
	
	* Tabulate results
	global esttabOptions "modelwidth(15) varwidth(30) star(* 0.1 ** 0.05 *** 0.01) nogaps noconstant not"	
	if "$sr_14to19" == "1" {
		* Start log
		cap log close sr_14to19
		log using "$path_logs/b1_sr_14to19.log", replace text name(sr_14to19)
		
		* Load estimations 
		estimates clear
		local estNumbers "14/19 1014/1019 10014 20014"
		local estimateFiles ""
		foreach e of numlist `estNumbers' {
			local estimateFiles`e'_u : dir "$root/savedEstimations/" files "`e'up_*.ster", respectcase
			local estimateFiles`e'_d : dir "$root/savedEstimations/" files "`e'down_*.ster", respectcase
			local estimateFiles`e'_a : dir "$root/savedEstimations/" files "`e'abs_*.ster", respectcase
			local estimateFiles`e'_f : dir "$root/savedEstimations/" files "`e'ff4_*.ster", respectcase 
			
			local estimateFiles`e' : list estimateFiles`e'_u | estimateFiles`e'_d
			local estimateFiles`e' : list estimateFiles`e' | estimateFiles`e'_a
			local estimateFiles`e' : list estimateFiles`e' | estimateFiles`e'_f
			local estimateFiles : list  estimateFiles  | estimateFiles`e'
		}
		
		foreach estimateFile of local estimateFiles {
			estimates use "$root/savedEstimations/`estimateFile'"
			local newEstName = subinstr("`estimateFile'", "_FPPI", "", .)
			local newEstName = subinstr("`newEstName'", ".ster", "", .)
			local newEstName = subinstr("`newEstName'", "LPM_FE", "", .)
			_eststo e`newEstName'
		}
		
		** Base
		capture noisily esttab e14up* e1014up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e14down* e1014down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ebpnew ${ageControl})
		//capture noisily esttab e14abs* e1014abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e14ff4_up* e1014ff4_up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e14ff4_down* e1014ff4_down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		//capture noisily esttab e14ff4_abs* e1014ff4_abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e15up* e1015up*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e15down* e1015down*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.ebpnew_residualP ebpnew ${ageControl})
		//capture noisily esttab e15abs* e1015abs*, ${esttabOptions} keep(L_x*_Px c.absGap_*#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e15ff4_up* e1015ff4_up*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e15ff4_down* e1015ff4_down*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		//capture noisily esttab e15ff4_abs* e1015ff4_abs*, ${esttabOptions} keep(L_x*_Px c.absGap_*#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		
		
		** TFE
		capture noisily esttab e16up* e1016up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ${ageControl})
		capture noisily esttab e16down* e1016down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ${ageControl})
		//capture noisily esttab e16abs* e1016abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.ebpnew_residualP ${ageControl})
		capture noisily esttab e17up* e1017up*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.ebpnew_residualP ${ageControl})
		capture noisily esttab e17down* e1017down*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.ebpnew_residualP ${ageControl})
		//capture noisily esttab e17abs* e1017abs*, ${esttabOptions} keep(L_x*_Px c.absGap_*#c.ebpnew_residualP ${ageControl})
		capture noisily esttab e16ff4_up* e1016ff4_up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm ${ageControl})
		capture noisily esttab e16ff4_down* e1016ff4_down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm ${ageControl})
		//capture noisily esttab e16ff4_abs* e1016ff4_abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.FF4_alt1_norm ${ageControl})
		capture noisily esttab e17ff4_up* e1017ff4_up*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.FF4_alt1_norm ${ageControl})
		capture noisily esttab e17ff4_down* e1017ff4_down*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.FF4_alt1_norm ${ageControl})
		//capture noisily esttab e17ff4_abs* e1017ff4_abs*, ${esttabOptions} keep(L_x*_Px c.absGap_*#c.FF4_alt1_norm ${ageControl})
		
		** Dummies
		capture noisily esttab e18up* e1018up*, ${esttabOptions} keep(*gapDummy* ebpnew* ${ageControl})
		capture noisily esttab e18down* e1018down*, ${esttabOptions} keep(*gapDummy* ebpnew* ${ageControl})
		//capture noisily esttab e18abs* e1018abs*, ${esttabOptions} keep(*gapDummy* ebpnew* ${ageControl})
		capture noisily esttab e18ff4_up* e1018ff4_up*, ${esttabOptions} keep(*gapDummy* FF4_alt1_norm ${ageControl})
		capture noisily esttab e18ff4_down* e1018ff4_down*, ${esttabOptions} keep(*gapDummy* FF4_alt1_norm ${ageControl})
		//capture noisily esttab e18ff4_abs* e1018ff4_abs*, ${esttabOptions} keep(*gapDummy* FF4_alt1_norm ${ageControl})
		capture noisily esttab e19up* e1019up*, ${esttabOptions} keep(*gapDummy* ebpnew* ${ageControl})
		capture noisily esttab e19down* e1019down*, ${esttabOptions} keep(*gapDummy* ebpnew* ${ageControl})
		//capture noisily esttab e19abs* e1019abs*, ${esttabOptions} keep(*gapDummy* ebpnew* ${ageControl})
		capture noisily esttab e19ff4_up* e1019ff4_up*, ${esttabOptions} keep(*gapDummy* FF4_alt1_norm ${ageControl})
		capture noisily esttab e19ff4_down* e1019ff4_down*, ${esttabOptions} keep(*gapDummy* FF4_alt1_norm ${ageControl})
		//capture noisily esttab e19ff4_abs* e1019ff4_abs*, ${esttabOptions} keep(*gapDummy* FF4_alt1_norm ${ageControl})
				
		
		* Close log
		log close sr_14to19
	}
	
	if "$sr_20to26" == "1" {
		* Start log
		cap log close sr_20to25
		log using "$path_logs/b1_sr_20to25.log", replace text name(sr_20to25)
		
		* Load estimations 
		estimates clear
		local estNumbers "20/25 1020/1025"
		local estimateFiles ""
		foreach e of numlist `estNumbers' {
			local estimateFiles`e'_u : dir "$root/savedEstimations/" files "`e'up_*.ster", respectcase
			local estimateFiles`e'_d : dir "$root/savedEstimations/" files "`e'down_*.ster", respectcase
			local estimateFiles`e'_a : dir "$root/savedEstimations/" files "`e'abs_*.ster", respectcase
			local estimateFiles`e'_f : dir "$root/savedEstimations/" files "`e'ff4_*.ster", respectcase 
			
			local estimateFiles`e' : list estimateFiles`e'_u | estimateFiles`e'_d
			local estimateFiles`e' : list estimateFiles`e' | estimateFiles`e'_a
			local estimateFiles`e' : list estimateFiles`e' | estimateFiles`e'_f
			local estimateFiles : list  estimateFiles  | estimateFiles`e'
		}
		
		foreach estimateFile of local estimateFiles {
			estimates use "$root/savedEstimations/`estimateFile'"
			local newEstName = subinstr("`estimateFile'", "_FPPI", "", .)
			local newEstName = subinstr("`newEstName'", ".ster", "", .)
			local newEstName = subinstr("`newEstName'", "LPM_FE", "", .)
			_eststo e`newEstName'
		}
		
		* US results
		** Pre 2007
		capture noisily esttab e20up* e1020up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e20down* e1020down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ebpnew ${ageControl})
		//capture noisily esttab e20abs* e1020abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e20ff4_up* e1020ff4_up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e20ff4_down* e1020ff4_down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		//capture noisily esttab e20ff4_abs* e1020ff4_abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e21up* e1021up*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e21down* e1021down*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.ebpnew_residualP ebpnew ${ageControl})
		//capture noisily esttab e21abs* e1021abs*, ${esttabOptions} keep(L_x*_Px c.absGap_*#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e21ff4_up* e1021ff4_up*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e21ff4_down* e1021ff4_down*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		//capture noisily esttab e21ff4_abs* e1021ff4_abs*, ${esttabOptions} keep(L_x*_Px c.absGap_*#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		
		** Low frequency
		capture noisily esttab e22up* e1022up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e22down* e1022down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ebpnew ${ageControl})
		//capture noisily esttab e22abs* e1022abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e22ff4_up* e1022ff4_up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e22ff4_down* e1022ff4_down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		//capture noisily esttab e22ff4_abs* e1022ff4_abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e23up* e1023up*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e23down* e1023down*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.ebpnew_residualP ebpnew ${ageControl})
		//capture noisily esttab e23abs* e1023abs*, ${esttabOptions} keep(L_x*_Px c.absGap_*#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e23ff4_up* e1023ff4_up*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e23ff4_down* e1023ff4_down*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		//capture noisily esttab e23ff4_abs* e1023ff4_abs*, ${esttabOptions} keep(L_x*_Px c.absGap_*#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		
		** High frequency
		capture noisily esttab e24up* e1024up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e24down* e1024down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ebpnew ${ageControl})
		//capture noisily esttab e24abs* e1024abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e24ff4_up* e1024ff4_up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e24ff4_down* e1024ff4_down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		//capture noisily esttab e24ff4_abs* e1024ff4_abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e25up* e1025up*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e25down* e1025down*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.ebpnew_residualP ebpnew ${ageControl})
		//capture noisily esttab e25abs* e1025abs*, ${esttabOptions} keep(L_x*_Px c.absGap_*#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e25ff4_up* e1025ff4_up*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e25ff4_down* e1025ff4_down*, ${esttabOptions} keep(L_x*_Px c.L_x*_Px#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		//capture noisily esttab e25ff4_abs* e1025ff4_abs*, ${esttabOptions} keep(L_x*_Px c.absGap_*#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
	
		* Close log
		log close sr_20to25
	}

		if "$sr_26" == "1" {
		* Start log
		cap log close sr_20to25
		log using "$path_logs/b1_sr_26.log", replace text name(sr_26)
		
		* Load estimations 
		estimates clear
		local estNumbers "26 1026"
		local estimateFiles ""
		foreach e of numlist `estNumbers' {
			local estimateFiles`e'_u : dir "$root/savedEstimations/" files "`e'up_*.ster", respectcase
			local estimateFiles`e'_d : dir "$root/savedEstimations/" files "`e'down_*.ster", respectcase
			local estimateFiles`e'_a : dir "$root/savedEstimations/" files "`e'abs_*.ster", respectcase
			local estimateFiles`e'_f : dir "$root/savedEstimations/" files "`e'ff4_*.ster", respectcase 
			
			local estimateFiles`e' : list estimateFiles`e'_u | estimateFiles`e'_d
			local estimateFiles`e' : list estimateFiles`e' | estimateFiles`e'_a
			local estimateFiles`e' : list estimateFiles`e' | estimateFiles`e'_f
			local estimateFiles : list  estimateFiles  | estimateFiles`e'
		}
		
		foreach estimateFile of local estimateFiles {
			estimates use "$root/savedEstimations/`estimateFile'"
			local newEstName = subinstr("`estimateFile'", "_FPPI", "", .)
			local newEstName = subinstr("`newEstName'", ".ster", "", .)
			local newEstName = subinstr("`newEstName'", "LPM_FE", "", .)
			_eststo e`newEstName'
		}
		
		** With weights 
		capture noisily esttab e26up* e1026up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e26down* e1026down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP ebpnew ${ageControl})
		//capture noisily esttab e26abs* e1026abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.ebpnew_residualP ebpnew ${ageControl})
		capture noisily esttab e26ff4_up* e1026ff4_up*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		capture noisily esttab e26ff4_down* e1026ff4_down*, ${esttabOptions} keep(L_x*_Pd c.L_x*_Pd#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		//capture noisily esttab e26ff4_abs* e1026ff4_abs*, ${esttabOptions} keep(L_x*_Pd c.absGap_*#c.FF4_alt1_norm FF4_alt1_norm ${ageControl})
		
		* Close log
		log close sr_26
	}

	//#changed20191202 (bootstrap)
	if "$sr_10014" == "1" {
		* Start log
		cap log close sr_10014
		log using "$path_logs/b1_sr_10014.log", replace text name(sr_10014)
		
		* Load estimations 
		estimates clear
		local estNumbers "10014 20014"
		local estimateFiles ""
		foreach e of numlist `estNumbers' {
			local estimateFiles`e'_u : dir "$root/savedEstimations/" files "`e'up_*.ster", respectcase
			local estimateFiles`e'_d : dir "$root/savedEstimations/" files "`e'down_*.ster", respectcase

			
			local estimateFiles`e' : list estimateFiles`e'_u | estimateFiles`e'_d
			local estimateFiles : list  estimateFiles  | estimateFiles`e'
		}
		
		foreach estimateFile of local estimateFiles {
			estimates use "$root/savedEstimations/`estimateFile'"
			local newEstName = subinstr("`estimateFile'", "_FPPI", "", .)
			local newEstName = subinstr("`newEstName'", ".ster", "", .)
			local newEstName = subinstr("`newEstName'", "LPM_FE", "", .)
			_eststo e`newEstName'
		}
		
		* Bootstrap test
		capture noisily esttab e10014up* e20014up*, modelwidth(15) varwidth(30) star(* 0.1 ** 0.05 *** 0.01) nogaps noconstant se keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP* ebpnew) rename(c.L_x_Pd#c.ebpnew_residualP_bs c.L_x_Pd#c.ebpnew_residualP) mtitles("Regular se" "Bootstrap se") title("Base, h12, up")
		esttab e10014down* e20014down*, modelwidth(15) varwidth(30) star(* 0.1 ** 0.05 *** 0.01) nogaps noconstant se keep(L_x*_Pd c.L_x*_Pd#c.ebpnew_residualP* ebpnew) rename(c.L_x_Pd#c.ebpnew_residualP_bs c.L_x_Pd#c.ebpnew_residualP) mtitles("Regular se" "Bootstrap se") title("Base, h12, down")
	
		* Close log
		log close sr_10014
	}




	
    * Close log
    log close summariseResults
}
